global KaiA KaiB KaiC m kxy0 kxyA Khalf k1 k2

% 1,2,3,4 represent T,D,S,U respectively

tspan=[0:0.5:200];
y0=[0.68;
    1.36;
    0.34;
    KaiC-0.68-1.36-0.34;
    0.00];                         % original concentration: T=0.68,D=1.36,S=0.34,U=1,SB=0 , according to the paper

[t,Y]=ode45('odefunction_Amax',tspan,y0);

KK=KaiB*ones(length(Y(:,5)),1);
K2=k2*ones(length(Y(:,3)),1);
R=(Y(:,5)./Y(:,3))./((k1*(KK-Y(:,5)))./(K2+k1.*Y(:,3)));

A=max(0,KaiA-2*m*Y(:,5));

% lamda is the ratio of k_TD to k_TU, lamda_p is the remained concentration of T
lamda=(kxy0(1,2)+kxyA(1,2)*A./(Khalf+A))./(kxy0(1,4)+kxyA(1,4)*A./(Khalf+A));
lamda_p=((kxy0(1,2)+kxyA(1,2).*A./(Khalf+A)).*Y(:,1))-((kxy0(2,1)+kxyA(2,1).*A./(Khalf+A)).*Y(:,2));
%figure;plot(t,lamda);legend('lamda');
figure;plot(t,lamda_p);legend('lamda^,')

N=24;                               %T
for j=1:4
n=j;                                %number of cycles
t0=15;                              %starting point
tx=1:1:n*N;                         %time point
lamda_p_int=zeros(1,N*n);
lamda_p_int(1)=1/12*(lamda_p(t0+1)+4*lamda_p(t0+2)+lamda_p(t0+3));              %Simpson's rule
for i=3:2:2*n*N-1
    lamda_p_int((i+1)/2)=lamda_p_int((i-1)/2)+0.5/6*(lamda_p(t0+i)+4*lamda_p(t0+i+1)+lamda_p(t0+i+2));   %Simpson's rule
end
figure;plot(tx,lamda_p_int);legend('lamda^,') 
hold on
end
hold off
